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Abstract. A novel procedure for georeferencing raster images of the Third 
Military Survey maps from 19th century (1876-1880) is proposed in the 
contribution. There were several attempts to design a proper transfor- 
mation model for georeferencing of these maps, but none of them have sat- 
isfactory resolved crucial problems of positional discrepancy between the 
old maps and contemporary maps (almost 150 m in reality). In the contri- 
bution, these problems are approached with the aid of three unique princi- 
ples. 



1 Huge number of control points was utilized to properly match the 
old maps to the contemporary coordinate system. The set of the control 
poi nts was val i dated by means of stati sti cal tests. 

2. A complex transformation procedure comprising of four partial se- 
quential transformation steps (rectification of paper shrinkage, reverse pro- 
jection of a map sheet onto ellipsoid, cartographic projection from the ellip- 
soid to plane, special elastic 2D transformation). The most important step is 
the special elastic transformation witch corrects inhomogeneous distortions 
of the old map sheets. Simultaneously, precision of an arbitrary point can 
be easilyesti mated. 

3. The transformation model is adjustable by a simple set of parame- 
ters with comprehensible meaning, namely standard deviation of the per- 
missible positional discrepancy. 



Due to these three principles the positional accuracy was reduced to a few 
meters (9 - 10 m) in reality, which corresponds to 0,4mm on the map sheet. 
Such a high accuracy is sufficient to produce a seamless map of the Czech 
Republic composed from the separate transformed map sheets. Therefore, 
any part of this composition can be compared to some up-to-date map 
source, eventually to another properly georeferenced old map. Accordingly, 
the proposed procedure was implemented as a web application on web 
server of the Research I nstitute of Geodesy, Topography, and Cartography. 
Consequently, any I nternet user can create overlays of a region of interest 
from several coverages and compare them to the content of the Third M ili- 
tary Survey maps, e.g. with the aid transparency. 

Keywords: georeferencing, accuracy analysis, elastic transformation, digi- 
tal image processing, historical cartography 

1. Introduction 

The motivation for georeferencing raster images of the Third Military Sur- 
vey maps stems from the ambition to make them available to the profes- 
sional public on Internetasawebmapservice(WMS). Broad availability of 
these old maps enables creating case studies of various branches, namely 
comparative history, countryside development or urbanism. These studies 
can greatly benefit from the ability to compare the content of the old maps 
with contemporary maps or other georeferenced old maps. An example of 
this comparison is overlaying of the maps with properly chosen transparen- 
cy. One of the most desired map for this kind of comparison is the Third 
Military Survey map in 1 : 25 000 scale (so called topographic section) be- 
cause of high level of detail. If the georeferencing is correctly designed with 
the aid of a suitable transformation procedure, it should be possible to ap- 
ply the procedure to the Special Map of the Third Military Survey in 1: 75 
000 scale as wel I . The Special M ap was bei ng used i n our cartographic prac- 
tice for over 70 year and many emissions of these maps were issued. For 
that reason these maps and the topographic sections are some of the most 
valuable cartographic sources depicting the development of the Czech coun- 
tryside and urban centres. 

Georeferencing the maps of the Third Military Survey of the Austro- 
Hungarian Monarchy to contemporary cartographic coordinate system is a 
long-term problem that has not yet been satisfactorily resolved. This prob- 
lem concerns primarily map sheets in 1 : 25 000 scale, where the positional 
differences according to contemporary maps are clearly visible due to their 
larger scale. Magnitude of these differences are so high that the maps' use 
in WMS is practically impossible. Common georeferencing techniques, e.g. 



matching only the corners of map sheets, achieve positional errors in range 
of 12 to 206 m in the extent of the Czech Republic. Such errors are even 
significantly greater than georeference errors of the the Second Military 
Survey, which is much older. These errors can be treated as a result of im- 
perfection of the then geodetic control and, simultaneously, shortcomings 
of the method of the Third M ilitary Survey maps' creation, namely convert- 
ing insular maps of Stable Cadastre into the topographic sections. A num- 
ber of papers deal with the analysis of these errors, e.g.(Cechurova 2009), 
(Cada 2006), (Cada 2005), (Krnoul 2010), (Molnar and Timar 2009), 
(Molnar and Timar 2011), (Seemann 2008), (Krnoul 2012). The sheer 
number of works attests to the i mportance of the problem. 



2. Positional differences between III. Military Survey 
and reality 

Figure 1 shows the positional differences between the maps of 1 1 1 . M ilitary 
Survey and reality. The arrows represent shifts on a selection of 4 246 trig- 
onometric points obtained by transformation with the use of map sheet 
corners as control poi nts. As we can see, the shifts from two disti net groups, 
each with a common global trend - one group covers the Bohemia and west- 
ern Moravia and the other covers the rest of Moravia. The average shift 
length is 108 m with standard deviation of 28 m. Most of the shifts, 2370 
i.e. 56%, fall into the 90- 130 m range. The maximum shift length is 206 m. 

Figure 2 shows the consequences of local positional differences when trans- 
forming the maps only with the help of map sheet corners. Contemporary 
National Map 1: 5 000 (deep black outlines of buildings) is overlapped with 
the transformed III. Military Survey map (gray and colour drawing). The 
mutual shift is so substantial that visual comparison of both overlapping 
maps i s al most i mpossi bl e. 

Because of these significant problems the best georeference precision 
achieved to this days has been in the range of several tens of meters (cca 40 
m) in reality (see (Cechurova 2009), (Cada 2006), (Krnoul 2010), (Molnar 
and Timar 2009), (Molnar and Timar 2011), (Seemann 2008), (Krnoul 
2012)). It isn't as much important why and how these errors emerged as 
finding a georeferencing method that would eliminate them as much as 
possible. 



to 




Figure! Positional errors on trigonometric points in the Third Military Survey 




Figure 2, Overlap of the III. M ilitary Survey map and contemporary map. Georef- 
erenci ng was done via transformation with the hel p of map sheets corners. 



3. The proposed method of georeferencing 

This paper proposes a new method for georeferencing the topographic sec- 
tions of the III. M ilitary Survey. This method takes advantage of unconven- 
tional approaches based on modern statistical methods for experimental 
data processing. These approaches build on three principles. 

1 It is necessary to use as many control points as possible for georef- 
erencing. The validity of the points must be verified by statistical 
tests. 

2. Paper distortion must be rectified and the original map projection 
must be respected. A special elastic transformation should be used 
to account for non- homogenous inaccuracy distribution in old maps 
without causi ng unacceptable deformations of the maps' content. 



3. The parameters of the transformation model should have clear 
meaning to enable fine-tuning the model to given situation. 

The proposed process consists of four steps. 

1 Elimination of paper distortion (correcting the sizes of map sheets 
to original non-distorted size). 

2. Projection of the map sheets on Bessel ellipsoid (a projection inverse 
to the original map projection) 

3. Map projection of Bessel ellipsoid into plane (projection into the 
current map projection of choice, e.g. Kfovak's projection). If the se- 
lected contemporary map projection is based on other ellipsoid than 
Bessel 's, it is necessary to transform the Bessel ellipsoid coordinates 
to the appropri ate el I i psoi d f i rst. 

4. Elastic transformation in plane (to correct the non- homogeneous 
map precision distribution). 

Each of these steps is an intermediate transformation. The final transfor- 
mation is a composition of all these intermediate transformations. The re- 
sult of appl i cation of the fi nal transformation on the source raster i mages of 
map sheets is digital images of the maps in the desired coordination system. 
From these images we can assemble seamless mosaic covering the whole 
area of the Czech Republ i c. 



4. Map distortion elimination 

Affine transformation is used to eliminate the map distortion, because the 
paper is usually shrunken differently in different directions. The decisive 
factor in this is the manufacturing process of paper that causes the distor- 
tion to have extreme values in approximately perpendicular directions. 

x = QX + r 

x ... input coordinates 

X ... output coordinates 

Q ... affine transformation matrix 

r ... coordi nate system ori gi n offset 

The corners of map sheets act as control points. Their coordinates were 
measured in source raster images. Their coordinates in target coordinate 
system can be calculated from Bessel ellipsoid parameters with the condi- 



tion of preserving the length of central meridian and lateral parallels of four 
topographic sections (see chapter 5). 

The affine transformation's equation for control points is: 



v t ... input coordinates residual of control point i 

The affine transformation matrix elements and coordinate system origin 
offset were estimated by the Least Squares method (Q,f ). With these pa- 
rameters known, the transformation formula for the map sheets' scans can 
be assembled as: 



5. Projecting the map sheets on Bessel ellipsoid 

The III. Military Survey uses piecewise Sanson- F I amsteed projection. In 
this projection, the surface of Bessel ellipsoid is divided into 30' x 15' (lon- 
gitude x latitude) sections and each section is projected into a plane by 
means of bilinear projection preserving the length of central meridian and 
lateral parallels of the section. This results in a trapezoid of appropriate 
size. The image of the 30' x 15' section is called section map sheet (in 1 : 75 
000 scale) and contains four topographic sections ( in 1 : 25 000 scale). 
Detailed information about Sanson projection can be found for example in 
(Boguszak 1961), (Cada 2006), (Seemann 2008), (Kriioul 2012). The area 
covered by a section sheet on Bessel ellipsoid is obtained my means of in- 
verse Sanson projection. In both directions the transformation can be de- 
fined with the help of so called plate. The source coordinates were normal- 
ized to (0,1) interval. 



x 



Xi+v t = QXi + r 

... input coordinates of control point i 
... resulting coordinates of control point i 



X = Q^ix-r) 
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F n (t) = l-n + (2n-l)t. 



v(To> T i) ■■■ general point of plate coordinates, T ,T t e (0,1) 

x k l ... coordinates of plate corners 

a it j ... pi ate edges -line segments between pi ate corners 



F n (t) --- hel per fundi on, n e {0,1}, t e (0,1) 

The meaning of symbols i,j,k,l is explained on Figure 3. 
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Figure 3. The meaning of indices of corners and edges of map sheets 

I ndices pairs in corners of the rectangles are indices k, I. The indices next to 
the arrows denote the pi ate edges, i.e. indices 

To solve the inverse transformation we only need to switch the source and 
target coordinates systems. However it is always necessary to normalize the 
source coordi nates to the (0,1) i nterval . 



6. Map projection of Bessel ellipsoid to plane 

We decided to use Kfovak's projection, because we wanted results in S- 
J TSK coordinate system. The projection equations can be found in a num- 
ber of publications and textbooks, e.g. (Kostelecky et al. 2010). The equa- 
tions were implemented in custom transformation program. 



7. Elastic planar transformation 

The elastic transformation in the S-J TSK plane was done with the use of 
collocation. This method allows to find a transformational relation between 
two coordi nate systems that takes i nto account both the control poi nt preci- 
sion and non- homogeneous distortion in the coordinate systems. The in- 
homogenity causes residual positional deviations on control points after 
Helmert transformation. The deviations cannot be justified by inaccurate 
measurement of control points' positions and therefore we can consider 
them a random variable describing the differences between elastic and 
similarity transformations. These differences cause nonlinear distortion of 
equidistant coordinate grid, creating the impression of elasticity (see Figure 



4). 



Figure 4. Illustrative example of equidistant coordinate grid distortion of 
a mapsheet caused by elastic transformation. The control poi nts are at the corners. 

These random deviations are statistically dependent. The exact dependency 
is characterized by covariance matrix of any set of points. The size of non- 
diagonal elements (so called covariances) diminishes with distance of the 
points in question. Both factors (control points precision and non- 
homogeneous distortion of coordinate systems) can be described by intui- 
tive statistical and geometric parameters: standard deviations of errors of 
control points position determination, standard deviation of difference be- 
tween elastic and similarity transformations in any point and a characteris- 
tic of distances between any two points in which their dependency mani- 
fests. The accuracy of the elastic transformation is estimated on the basis of 
these parameters. Control points are trigonometric points with known S- 
J TSK coordinates depicted on the map sheets. 

The col location method is based on the following simultaneous equations: 



w' ... target coord i nates 

p, q ... similarity transformation coefficients 

(p,(Pi ... random differences between elastic and similarity transfor- 
mations 



w' = p + qw + cp 

w'i + s'i=p + q(w t + Ei)(pi 



w 



... source coordinates 




w 



... source coordinates of control point i 

... target coordinates of control point i 

... positional errors of coordinates of control point i 



w',p,q ... unknown variables 



All of these variables have complex values. Complex arithmetic is also used 
to calculate the values of unknown variables. This (complex) approach sig- 
nificantly simplifies the computer solution of the simultaneous equations. 
To solve the inverse transformation we only need to switch the source and 
target coordination systems. 

8. The final composite transformation 

The final transformation is a composite projection consisting of sequential 
application of the individual transformations. It is a fairly complex trans- 
formation and a demanding calculation which makes it unfeasible to calcu- 
late exact values for every pixel. Because of that suitable grid of node points 
ischosen. Pixel values are computed according the exact equations for these 
node points. The values for the rest of pixel are obtained via bilinear inter- 
polation. This approach allows to substantially reduce the computing time 
while preserving the required pixel precision. 

The result of the final transformation is a digital image of the source map in 
contemporary map projection coordinate system (here it is S-J TSK). These 
images can be arranged into a seamless mosaic covering the whole Czech 
Republic. An exampleof this mosaic ison Figure5. 




Figure 5. Topographic sections mosaic covering the northern part of the Czech 
Republic. Blue dries represent control points. 



The original scans of map sheets were transformed by nearest neighbour 
method which associates pixel in resulting images with the colour of its 
nearest neighbour in source images. However an inverse transformation, 
i.e. a sequence of inverse component transformations arranged in reverse 
order, is needed for this. It was therefore necessary to implement two ver- 
sions of the transformation: forward and backward. Forward transfor- 
mation was used to estimate the accuracy of resulting image, backward 
transformation was used to determine the pixel colour in resulting image 
(with the nearest neighbour method). An alternative to the nearest neigh- 
bour is bilinear or bicubic interpolation. These methods provide smoother 
and more visually pleasing resulting image, but they do not contribute to 
better accuracy because they only manifest on a sub-pixel level (0.06 mm). 
To improve the visual quality of the resulting image one of these method 
will be likely selected in the future based on detailed comparisons. 




Figure 6. An example of seamless connection of 4 map sheets. The corners' con- 
nection is in the middle of the red ellipse. 




Figure 7 - Overlap of the 1 1 1 . M i I itary Survey map and contemporary map. Georef- 
erenci ng was done by the proposed method. 



The results of the final transformation can be seen on Figures 6 and 7. 
Figure 6 shows the connection of four neighbouring map sheets (3853/2, 
3852/ X 3853/ 4, 3852/4). The connection is in the center of a red ellipse. 
Control points are represented by blue circles and red numbers. The first 
five digits designate the maps sheet's signature. 

Figure 7 shows and overlay of contemporary vector map data (buildings 
layer of National Map l : 5 000, provided by CUZK) over the III. M i I itary 
Survey map sheet. The contemporary map layer is shown as deep black 
buildings' outlines over the colour background of the III. Military Survey 
map. Unlike the situation on Figure 2, the two maps here correspond to 
each other very wel I . 



9. The method for data acquisition 



The topographic sections of the III. M ilitary Survey in 1 : 25 000 scale were 
acquired from several sources to cover as much of the Czech Republic's area 
as possi ble. Some of them were found i n the map col lection of The I nstitute 
of H i story, Academy of Sciences of the Czech Republ ic, v.v.i . and scanned i n 
VUGTK, v.v.i. Some of the missing colour originals had to be replaced by 
black-and-white copies. Out of the 376 topographic sections there is 234 in 
colour, 133 map sheets are black-and-white and 9 map sheets are still miss- 
ing. We are still looking for the missing map sheets in various map collec- 
tions i ncl udi ng private col lection. The current coverage of the Czech Repub- 
liccan be seen on Figure8. 




Figure 8. Overview of scanned topographic sections 



To georeference the maps it is necessary to have some control points and 
obtain images coordinates of map sheets corners. Trigonometric points and 
churches with known coordinates were selected as control points. Details 
about the number of various objects' representation in the control points' 
set is in Table! 



CP type 


CP description 


Simplified map 
symbol 


Real map symbol (on 
map) 








Colour 


Black- 
and- 
white 


TB 


Trigonometric point 


A 






Kl 


Church - type 1 


i 










K2 


Church - type 2 












K3 


Church - type 3 


t 


i 




K4 


Church - type 4 


i 


% 





Table 1 Control point types 



As many control points as possible were used for each topographic section. 
Nine sections do not contain any usable control point. These sections depict 
border areas of the Czech Republic and show only a small part of the Czech 
territory. The rest of map sheets contain between 1 and 28 control points 
each. The image coordinates of control points were measured in geodetic 
program Kokes. The identification of control points was often very difficult 
because of the quality of the map print especially on the black-and-white 
maps. On some topographic sections even positions of the section corners 
could not be reliably determi ned. There were three cases of map sheets with 
the corners cut away, probably as a result of mistakes and unskilled manip- 
ulation with the paper originals. Fifteen sections were reduced to half- or 
quarter the usual size, likely to save paper when drawing small areas. And 
finally 6 pairs of topographic sections were joined into a single, larger sec- 
tion. 

The joined topographic sections: 3654/2 +3554/4 

3850/1 + 3750/3 
3857/1 + 3857/3 
4455/2+4455/4 
4456/4+4556/2 



4457/ 3 + 4557/ 1 



The whole process of finding the control points and measuring their 
coodi nates was done manually and the result was 4526 control points and 
about 1400 corners points. 




Figure 9. Topographic sections without control points - red rectangles 
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Table 2. Missing corners and irregular size of topographic sections 




Figure ID. Positional errors on check points 



10. Conclusion 



The proposed method of georeferencing contributes significantly to the us- 
ability of the III. M ilitary Survey maps for practical applications ( such as 
comparative history, countryside development urbanism, planning, tour- 
ism, etc.). This contribution is based on reduction of errors of the old map 
placement into contemporary coordinate system. It is shown, that the posi- 
tional discrepancy of map elements that did not change over time can be 
reduced so much as to make visual comparison possible. 

Trigonometric points and church towers were used as control points. The 
positional accuracy was evaluated with the use of other 958 check points. 
The positional errors on these points are shown on Figure 10. The resulting 
positional deviation for the check points is 9.1 m. This is a substantial im- 
provement compared to previous gereferencing attempts. This is mainly 
because of the step 4 of the georeferenci ng process, i .e. the elastic transfor- 
mation. The elastic transformation in conjunction with a large quantity of 
control points allows seamlessly join the individual map sheets, correct the 
local inaccuracies of the map (deformation and distortion) and properly 
estimate the final positional accuracy of the transformation. The formal 
correctness of the accuracy estimation is ensured by the statistical proper- 
ties of collocation. 

The purpose of georeferenci ng the 1 1 1 . M i I itary Survey was achieved and the 
results can be published as WMS on server mapy.vugtk.cz. Internet users 
will be given the ability to use this service in their own applications and 
studies concerni ng the study of content of the III. Mi I itary Survey maps and 
its comparison with current reality. 

This paper was created as part of project n. DF1P01OW021 "Cartographic 
sources as a cultural heritage. Research of the new methods and technolo- 
gies of digitalisation to enable access and use of the old maps, plans, atlases 
and globes" under the auspices of the Program for applied research and 
development of national and cultural identity. 
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